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Formation flying is an enabling technology for many future space missions. This pa- 
per presents extensions to the equations of relative motion expressed in Keplerian orbital 
elements, including new initialization techniques for general formation configurations. A 
new linear time-varying form of the equations of relative motion is developed from Gauss’ 
Variational Equations and used in a model predictive controller. The linearizing assump- 
tions for these equations are shown to be consistent with typical formation flying scenarios. 
Several linear, convex initialization techniques are presented, as well as a general, decen- 
tralized method for coordinating a tetrahedral formation using differential orbital elements. 
Control methods are validated using a commercial numerical propagator. 

I. Introduction 

I N spacecraft formation flying mission design, it is generally true that the control of relative spacecraft 
states is more important than the control of absolute states. In addition, knowledge of the relative states 
of spacecraft in a formation is often far more accurate than knowledge of the formation’s absolute state. 
For these reasons, formation control objectives are typically focused on controlling the relative states of the 
spacecraft. Variants of Hill’s and Lawden’s equations of relative motion have been used for online planning 
and control.^ Both of these approaches linearize the nonlinear dynamics of orbital motion about a reference 
orbit. However, the linearization approach taken in those cases restricts the separation distances of the 
satellites in the formation. For larger separation distances, the equations of motion can no longer be used 
to cancel relative drift rates (initialization) or to accurately predict the effect of inputs (control) In the 
case of the MMS mission, four spacecraft will be placed in a tetrahedron-shaped relative configuration that 
will have sides ranging between 10-1000 km. These separations far exceed the distances for which Hill’s and 
Lawden’s models are valid, even with the correction terms introduced in Ref. [24] 

Another approach that often appears in the literature is formation control using Gauss’ Variational 
Equations (GVEs)ra[B GVEs have been used for many years to account for perturbations in orbits 
arising from drag and Earth oblateness effects, as well as to design Lyapunov and fixed impulse control 
systems.^ GVEs are convenient for specifying and controlling widely separated formations because they 
are linearized about orbital elements, which are expressed in a curvilinear frame in which large rectilinear 
distances can be captured by small element perturbations. In addition, the GVEs provide a computationally 
simple way (no frame rotations are required) to obtain linearized dynamics about the orbits of each spacecraft 
in the formation. This bypasses the linearization error created by representing the entire formation in a single 
rectilinear frame, which was the approach used in Ref. [5j The use of GVE dynamics as opposed to Hill’s 
dynamics incurs the cost of computation associated with the use of multiple sets of time-varying equations 
of motion. Specifying a formation’s relative geometry in terms of differential orbital elements is an exact 
approach that does not degrade for large spacecraft separations. However, the advantage of using GVEs for 
control could be reproduced by using a separate Lawden’s frame for each spacecraft in the formation while 
still using orbital element differences to represent the formation relative geometry. Given that a nonlinear 
transformation and rotation is required to switch between a Lawden’s frame and orbital element differences, 
and that GVEs are already linearized in an orbital element frame, it is both simpler and computationally 
more efficient to use orbital element differences to specify the formation configuration and GVEs for control. 
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Several research groups have proposed control laws for formation-flying spacecraft that use GVEs to 


design impulsive thrusting maneuvers for orbit correction. For example, Ref. 14 proposes an general orbit 
correction scheme that uses GVEs to develop four impulsive thrusts that are applied at fixed points in an 
orbit. This four-impulse method is not guaranteed to be fuel-optimal and the approach presented in this 
paper consistently produces trajectories that require less fuel to accomplish identical goals. A method of 
producing optimized four-impulse plans for very-low eccentricity orbits is presented in Ref. |10[ but this 
approach does not extend to the higher eccentricities required for MMS missions. Another method based on 
GVEiP^ allows optimized planning for low Earth orbits, but only permits optimization over a single impulsive 
thrust, guaranteeing that the solution will be sub-optimal in many cases. In addition, this approach is 
only derived for correcting errors in semimajor axis, eccentricity, and inclination. Another approach to using 
GVEs for formation control is to derive a continuous proportional-derivative controller satisfying a Lyapunov 
equation 000113 

This paper presents a formation flying spacecraft control approach using GVEs as the dynamics in a 
model predictive control system. A novel aspect of this approach is the use of the GVEs to optimize the 
control inputs, or equivalently the spacecraft trajectory, over a set time-horizon. Plans are regularly re- 
optimized, forming a closed-loop system. By extending previous planning approaches to use GVEs, we can 
optimize the plans for spacecraft in widely-separated, highly elliptic orbits. Results are presented to show 
that the GVE-based planning system is more fuel-efficient than both the four-impulse method in Ref. [14] 
and the single-impulse optimized method in Ref. 8 In addition control optimized online has the advantage 
of being capable of handling many types of constraints, such as limited thrust capability, sensor noise 
robustness, and error box maintenance.^ This paper also presents a method for optimizing the initialization 
of the tetrahedron formation using relative orbital elements, which allows for rotation and translation of the 


formation. We also extend the virtual center approach to formation flying in Ref. 11 to GVEs and present 
a decentralized implementation. 

Gauss’ Variational Equations (GVEs) are derived in Ref. [16 and are reproduced here for reference 
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where the state vector elements are a (semimajor axis), e (eccentricity), i (inclination), (right ascension of 
the ascending node), u> (argument of periapse), and Mq (mean motion). The other terms in the variational 
expression are p (semi-latus rectum), b (semiminor axis), h (angular momentum), 6 (argument of latitude), 
r (magnitude of radius vector), and n (mean motion). All units are in radians, except for semimajor axis 
and radius (meters), angular momentum (kilogram • meters 2 per second), mean motion (1/seconds), and 
eccentricity (dimensionless). The input acceleration components u r , ug, and Uh are in the radial, in-track, 
and cross-track directions, respectively, of an LVLH frame centered on the satellite and have units of meters 
per second 2 . The form of the GVEs can be more compactly expressed as 

e = A(e) + B(e) u (2) 

where e is the state vector in Eq. [lj B(e) is the input effect matrix, u is the vector of thrust inputs in 
the radial, in-track, and cross-track directions, and A(e) = ( 0 0 0 0 0 \/n/a 3 ) T , where p is the 

gravitational parameter. 


II. Previous Approaches to Control Using GVEs 

A common approach when basing control on GVE dynamics is to use a nonlinear PD-based Lyapunov 
regulator of the form u _ ( 3 ) 

where £ is the current orbital element offset from the reference orbit e and K is a constant positive definite 
matrix 00EE3 Control algorithms of this type have been shown to be asymptotically stable in most casesj^ 
but belong to a class of control systems that fire continuously. Continuous firing is generally not desirable 
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for space missions because it is often disruptive to the science mission, it typically must be coupled with 
attitude maneuvers, and it expends fuel (nonreplenishable aboard a spacecraft) continuously. 

Another approach to differential element control is presented in Ref. |14| In that approach, it is observed 
that the GVEs decouple at several points during an orbit. By exploiting the decoupling points, an algorithm 
requiring a maximum of four impulsive thrusts is proposed. This approach is simple in an algorithmic sense, 
but requires a fixed time period ( z.e one orbit) to correct state errors and is not guaranteed to be fuel- 
optimal (or even near- fuel-optimal) . Ref. [8] presents another method of formation control based on GVEs 
that uses a single corrective thrust that is optimized nonlinearly. Although this method is guaranteed to find 
the optimal single-thrust correction for an arbitrary time period, it is not guaranteed (or likely) to find the 
optimal multiple-thrust correction. In addition, this approach is restricted to use in low Earth orbits and 
is only designed to correct errors in semimajor axis, eccentricity, and inclination. An approach presented in 
Ref. [15] uses a pseudo-inverse to the GVE control effect matrix to calculate a single corrective impulse. This 
approach is not guaranteed to be fuel-optimal for any cases and is not accurate for correcting position errors. 

In contrast, this section describes a control law that generally does not fire continuously and makes 
explicit its objective to minimize fuel use. The control approach utilizes the linearized relative dynamics of 
Gauss’ Variation Equations to optimize the effects of arbitrarily many inputs over a chosen planning horizon. 
The next section examines the validity of using a linearized relative form of the equations of satellite motion 
based on GVEs. 


III. Relative Orbital Elements and Linearization Validity 


In a formation, the orbital element state of the zth satellite is denoted e^. The states of the vehicles in 
the formation can be specified by relative orbital elements by subtracting the state of an arbitrarily chosen 
spacecraft in the formation (ei) 

Sei = e.i - ei (4) 

For a desired orbit geometry, a set of desired relative elements, Sedi will specify the desired state e* of each 
spacecraft in the formation 

&di - t (b) 

The state error for the zth spacecraft in the formation, Ci, is then defined as 

Ci ~ ^i ® di — &&di (6) 

The form of Gauss’ Variational Equations in Eq. [T] is for perturbations of orbital elements. To reformulate 
these equations for perturbations of relative orbital elements,^ the GVEs for and edi are placed together 

Ci ^i ® di 4l(e.j) A{(5di) 4“ B (e.i)lli (7) 


where the term B(edi)vLdi has been excluded because thrusting only occurs at the point of the spacecraft, 
not at its desired state, which is assumed to be a constant Keplerian orbit. The unforced dynamics can be 
linearized by introducing the first order approximation 


A( e i) — A(edi) 


dA . dA 

„ . (Gl “ Gdl) ~ 


Ci = A*(e d i)Ci 


( 8 ) 


where the matrix A* (edi) is all zeros except for the lower-leftmost element, which is — 3rz/2a. With this 
approximation, the differential GVE expression can be rewritten as 


Ci — A* (edi)Ci + B(ei)\ii — A*(edi)Ci + B(edi + Ci) u * 


( 9 ) 


In this case the control of the relative error state, Ci , is nonlinear, because the control effect matrix B is a 
function of the state. Ref. [13] accounts for this nonlinearity in a continuous nonlinear control law that was 
proven to be asymptotically stable. The control approach developed in this section uses linearized dynamics 
to predict the effect of future control inputs. Linearizing the matrix B in Eq. [9] yields 


6 


A*(edi)Ci + ( B(edi) + 


Ci 


A*(edi)Ci + B(edi)ui + [B* (e^jj^u,; 


( 10 ) 
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where the term B*(e d i) is a third rank tensor and the quantity B*(e d i)Ci is a matrix with the same dimensions 
as B(edi). For convenience, define 

AB = B*(e di )Ci (11) 


resulting in the new state equation 


Ci — ^*( e dz)Ci + ( B(edi ) + AS) Uj 


( 12 ) 


Note that if A B is much smaller than B(e d i), then the first order term can safely be ignored, yielding the 

linearized system ■ , 

C i = A ( Gdi)Ci + H( e di) u i (13) 

which can be controlled using any one of a variety of linear control techniques, including the model predictive 
controller discussed in Section UYI 

The critical requirement for linear control and planning is that the term A B holds much less influence 
on the state than the term B{e d i). However, A B is a linear function of the state error £ * , which can be 
arbitrarily large. The amount of acceptable error due to linearization will be a function of the mission 
scenario, but the linearization assumption will typically only be valid for small values of the state error. A 
bound on the magnitude of this error can be developed by comparing the induced norm of the difference 
between the control influence matrix at its desired state, B{edi) and at the actual position of the spacecraft, 
B(ei). In Eq. 11 the first order approximation of this term was defined as A B. In the following examples, 
A-Bt rue , which is defined as 

A B true = B{ei) - B(e di ) (14) 

and will be calculated numerically. The cutoff point of acceptable linearization error is when the norm of 
A B exceeds some (mission dependent) fraction of the norm of B{e d i). 

To investigate this cutoff point, the following examples consider many random values of Q in the set 
1 1 Ci 1 1 2 = t and calculate AB true . The A B true with the largest 2-norm will be used to test the validity of the 
linearization for a given r. This procedure is repeated for multiple r to find the largest HCII 2 for which the 
linearization is considered valid. 


Example: Low Earth Orbit — An example low Earth orbit is 

e di = ( 1.08182072 0.005000000 0.610865238 2t r 7 r 3.82376588 


(15) 


where the first element, semimajor axis, is normalized by the Earth’s radius, making the orbital element 
vector dimensionless. The matrix corresponding B{e d i) is 

( -5.6794478 1808.6011 0 \ 


B(e di ) = 


-0.000082308780 -0.00020502572 


0 


0 

0 

0.020528419 


0.00010304404 

0.00014406293 


-0.032987976 -0.00011800944 


\ -0.020792326 0.032987564 


(16) 




where || ^(e^i) || 2 = 1808.61. The effect of perturbing e d i for a given norm bound on Q is shown in Figure]!] 
The figure shows that for a linearization validity cutoff of 0.1, where || AB(e di , C)true||2 < 0.1 1| i3(e di ) true || 2 , 
can be achieved by observing the bound ||£j |2 < 0.08. This choice of bound allows for orbital element 
perturbations that equate to rectilinear distances on the order of 275 kilometers and velocities on the order 
of 416 meters per second. Since typical error box sizes for LEO formation flying missions are between 10-1000 
meters in size, the linearization should be a valid approximation for most missions of interest. 

Example: Highly Elliptical Earth Orbit — One motivation for using GVEs as the linearized dynamics 
in a planner is recent interest in widely spaced, highly elliptical orbits An orbit of this type is 


e di = 6.59989032 0.818181000 0.174532925 2ir 0 1 r 


( 17 ) 
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Fig. 1 : Effect of Orbital Element Perturbations on 
the ABtrue Matrix for a LEO Orbit 



Fig. 2 : Effect of Orbital Element Perturbations on 
the A£?true Matrix for a HEO Orbit 
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(18) 


Repeating the same procedure used for the LEO case, it is determined from Figure [2] that for a 10% 
linearization validity cutoff, the bound which must be observed is ||C ||2 < 0.034. In this case, the bound 
on II^H 2 corresponds to rectilinear distances of approximately 460 kilometers and velocities of 18.5 meters 
per second. As in the LEO case, these distances are far larger than expected error box sizes (i.e., expected 
distances over which control would be planned). Unlike the LEO case, error boxes for widely-separated 
missions, such as MMS, may be much larger than 10 meters to a side, even approaching 10’s of kilometers. 
The 10% cutoff ensures that error boxes of up to half the baseline are acceptable. In an actual mission, it is 
unlikely that error boxes would be large enough to impact the linearization validity, because it would allow 
the regularity of the tetrahedron geometry {i.e., the science data quality) to be significantly diminished. 


IV. Model Predictive Control Using GVEs 

Reference [5] showed that given a valid set of linearized dynamics and a desired trajectory, a model 
predictive controller for a spacecraft formation can be designed that allows for arbitrarily many convex 
terminal and intermediate state conditions, as well as sensor noise robustness requirements. This controller 
is implemented on each spacecraft in the formation and it is using a linear programming formulation. The 
general form of the optimization performed by the controller is 

min||u||i subject to Au < b (19) 

where u is a vector of potential control inputs and the matrix A and the vector b are formed based on the 
input dynamics and problem constraints. 

In order to use the linearized GVE-based dynamics developed in Eq. [13] in the MPC formulation, the 
dynamics are discretized using a zero order hold assumption according to the procedure described in Ref. |4] 
yielding the discrete form 

Ci(k + 1) = A d (e d i)£i(k ;) + B d (e d i)\ii(k) (20) 

where k is the current time step, A d (e d i) is the system’s parameter-varying state transition matrix, and 
Bdi^di) is the discretized form of the linearized GVE control input matrix B(e d i ). 
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V. Comparison to Another GVE-based Impulsive Control Scheme 


Section [II] describes a number of GVE-based control schemes. The controllers based on impulsive ap- 
proaches are those presented in Ref. [8! and Ref. 14 The single-thrust controller in Ref. |8]is restricted to use in 
LEO and can only control semimajor axis, eccentricity, and inclination. Ref. [14] describes a controller which 
uses four impulses over the course of an orbit to correct arbitrary orbital element perturbations. Because of 
its more general applicability, this section will compare the 4-impulse approach to the MPC-based approach 
presented in Section |IV| Both schemes are designed to drive the elements of a state error Q to zero over a 
fixed time interval. The four-impulse approach has not been presented in the context of performance criteria 
( e.g trajectory or terminal error boxes, robustness to disturbances) or constraints ( e.g ., maximum thrust 
level), so the comparisons in this section will use an MPC controller formulation that minimizes fuel use 
while driving the error state to zero in a fixed time and has no other constraints. 

The 4-impulse approach developed in Ref. [I4]can be summarized in four steps to be taken over the course 
of an orbit. When the angle of latitude, 9 , is 0 or n radians, implement a velocity change (impulsive thrust), 
A Vhi in the cross-track direction of an LVLH frame centered on the spacecraft to cancel the inclination error 
component of Q 

h a- 

A y hi = vAi 

r cos 9 

When the angle of latitude, 6 , is 7r/2 radians, implement a velocity change, A Vh n in the cross-track direction 
to cancel the ascending node error 

. h sin i — ^ 

Av hn = 

rsinO 


At perigee and apogee, implement Au 
of perigee and mean anomaly errors 

'(1 + e) 2 


and Av ra , respectively, in the radial direction to cancel the argument 


na 

Av r = 

r P 4 


V 


(Aw + Af 1 cos i) + AM 


na 

Al ’r« = ~4 


P — (Aw + AH cos i) + A M^j 


Also at perigee implement A vg p and at apogee implement A vg a in the in-track direction, to cancel the 
semimajor axis and eccentricity errors 




naij f A a 
~4~ l ”o" 



A Vg a 


nag f A a Ae \ 

~ 4 ~ V“a~ ~~ 1-eJ 


Using the notation and the HEO reference orbit from Section |III[ the following example compares the 
MPC method with the control approach reviewed in this section. For the state error 

Ci = ( 1CT 9 1CT 7 10 -7 1CT 7 1CT 7 10" 7 ) T (21) 

the 4-impulse method requires 1.42 mm/s of fuel to correct the state error over the course of an orbit and 
the MPC method requires 0.549 mm/s of fuel. For this example, the model predictive controller was given 
a full orbit time horizon. However, the same control objective could have been achieved in less time, but 
using more fuel. 

A series of 1000 orbital element state error vectors, Q, were generated in which each perturbed element 
was a random number between ±10 -6 . For each of the error vectors, both control methods were used to 
generate plans for eliminating the error. On average, the MPC maneuvers only required 51% of the fuel 
required by the 4-impulse maneuver. 


VI. General Drift-free Tetrahedron Initial Conditions 

The MMS mission will require spacecraft to form very large tetrahedron shapes, while simultaneously not 
drifting with respect to one another. Drift-free designs based on Hill’s and Lawden’s equations are valid only 
for formations with short baselines, because of the linearization assumptions inherent in the specification 
of the frames in the derivation of the dynamics. For any group of spacecraft, a no-drift requirement is 
equivalent to requiring that all spacecraft have the same orbital energy, which is also equivalent to stating 
that all spacecraft have the same semimajor axis. For a formation specified in differential orbital elements, 
this is the same as requiring that the desired differential semimajor axes for all spacecraft in the formation be 
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zero. Thus, in differential orbital elements it is trivial to design a drift-free formation, however, the curvilinear 
nature of the elements makes describing and manipulating general tetrahedron shapes complicated. 

One approach is to begin with tetrahedron coordinates in a rectilinear frame (such as LVLH or ECEF) 
and then convert those coordinates into orbital element perturbations. However, producing the desired 
orbital element differences for a tetrahedron requires knowledge of the full relative state in a rectilinear 
frame, including both position and velocity. If the only constraint on the formation geometry is that the 
satellite positions form a tetrahedron shape, then the velocity must be selected based on additional criteria. 

A regular tetrahedron with sides of length s can be described in a rectilinear frame using the coordinates 
(given in x, y , z triples) 


Sati = (0,0, fs) 


Sat 2 = (0, ^s,--^s) 


Sat 3 = (0,-^s, — ^-s) 


Sat 4 = (^s,0,0) 


If the rectilinear frame is LVLH, the coordinates can be transformed into differential orbital elements using 

the first order approximation . „ , . , . . 

Se di = M( ei )x dl (22) 


where e 4 is the reference orbit that the differential elements of the vector 5e di are described with respect to, 
and x.di is a relative state vector in the LVLH frame centered on the absolute orbit ei. The elements of the 
6x6 rotation matrix M{e i) are known analytically and can be found in Appendix G of Ref. 13 It should be 
noted that these projections from LVLH coordinates to differential orbital elements can be used to establish 
three and six dimensional “box” constraints on spacecraft position and velocity of the type discussed in 
Refs. [5l and [MI 

The desired LVLH state of satellite i can be expressed as the concatenation of a position vector p, and 
a velocity vector v,, each with states in the x, y , and z directions. 



(23) 


To find velocity vectors that satisfy the no-drift requirement, the system 


Se di = M{e l ) 



V i = 1 . . . n 


(24) 


must be solved with the additional constraint that the semimajor axis element of 5e d i be equal to zero for 
all spacecraft in the formation (n = 4 for a tetrahedron formation). In this system, the elements of Se d i and 
Vi are allowed to vary, while the matrix M and the vectors x,; are determined by the reference orbit and the 
tetrahedron geometry, respectively. The system has 9 n variables and 7 n constraints and will, therefore, have 
many possible solutions. In Ref. [18] the velocity magnitude is chosen to achieve the no-drift condition and 
the ECI velocity direction of each spacecraft in the formation is chosen to match the direction of reference 
orbit’s velocity vector. This method will succeed in creating a drift-free formation, however it does not take 
into account the states of the spacecraft in the formation immediately prior to initialization. 

The approach presented here minimizes the size of the maneuvers that would be required to create the 
desired formation. The approach selects initial conditions, Se di , that are closest to the current differential 
states, (ie,;, of the spacecraft in the formation and that will minimize the state error, Q, across the formation 
at the start of the initialization. For the entire formation, this criterion becomes 


min VllWi^-fe,)!!! (25) 

OGdi^di Vi=l...nr — : 

i=l 

where \V t are weighting matrices that represent the expected fuel-cost of changing orbital elements (obtain- 
able from the GVEs). Allowing different Wj for each spacecraft enables the formation design to take into 
account factors such as fuel-weighting to extend overall mission duration, similar to the approach used in 
choosing the virtual center in Section ESI The use of a 1-norm is appropriate in this case, because the 
distance that a given element must be changed is the absolute value of the difference between that element’s 
current and desired state. This approach is similar to the optimization used in Ref. [3, in which drift-free, 
minimum maneuver constants of integration were found for Lawden’s Equations. Next, the optimization is 
expanded by exploiting specific aspects of the MMS mission science goals. 

The quality of the shape of the regular tetrahedron^] largely determines the value of the science data 
recovered by a mission such as MMS. In choosing the initial conditions for a regular tetrahedron-shaped 


a Tetrahedron quality is discussed in Ref. 1 24 
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formation, certain quantities that do not affect shape quality such as the scale, position, and orientation of 


the tetrahedron can be considered degrees of freedom in the optimization described by Eqs. 24 and 25 By 
optimizing the additional degrees of freedom, tetrahedron-shaped initial conditions can be found that require 
smaller maneuvers to achieve from the current formation state. Scaling the tetrahedron shape equally in 
three dimensions introduces a single scalar variable s and allowing translation in each orthogonal direction 


of the LVLH frame introduces the variables t x , t y , and t z . The new constraint set is 


Se di = M(e i) 



(lOOOOo) 6e di = 0 V i = l...n 


(26) 


where t = ( t x t y t z 0 0 0 ) T and the second constraint forces a no-drift condition by ensuring that 

the relative semimajor axes, the first element of each Se di vector, are zero. These constraints can be combined 
with the objective in Eq. [25] and formulated as a linear program. In addition to the geometric and no-drift 
constraints, limits on the desired differential angle state variables are required to ensure that they remain 
within ± 7 r and on eccentricity and semimajor axis to ensure that the spacecraft remain in Earth orbit. 

For elliptical orbits, it is not possible to maintain constant relative geometry between satellites for all 
points in the orbit. Instead, a single position in the orbit must be chosen for the satellites to form a 
tetrahedron. The mean anomaly at the time of the tetrahedron geometry will be M t . When formulating 
the optimization in Eqs. |25| a nd |26| the current differential element vectors must be propagated forward 
using the A* matrix in Eq. [9] to the mean anomaly M t . In addition, the reference orbit used to compute the 
matrix M(e i) should have a mean anomaly set to M t . 

There are several limitations to this optimization approach. First, it does not optimally assign spacecraft 
to positions in the tetrahedron (a formulation that does this is possible using mixed integer linear program- 
ming or network LF023). Second, the optimization posed here does not optimally orient the tetrahedron in 
three space. Introducing a rotation matrix dependent on three Euler parameters would create a nonlinear 
optimization. This limitation could be bypassed by creating a spherical lattice about the orbit and opti- 
mizing the formation rotated once for each point on the lattice. After performing all the optimizations, the 
desired state corresponding to the rotation with the lowest cost would be chosen. Although this approach 
requires a preset number of optimizations (possibly many depending upon the degree of rotation resolution 
desired), the optimizations are small linear programs, which complete in a fraction of a second. 

An alternative form of Eq. [26] can be written to allow for small rotations using the linearized form of a 
three-dimensional rotation matrix 17 

1 e z -Oy 

1 Or I x = Rx (27) 


On 


-Ox 


1 


where 0 X is the rotation about the x axis, 0 y is the rotation about the y axis, 0 Z is the rotation about the z 


axis, x is an arbitrary LVLH position vector, and x' is the vector x after having been rotated. Using Eq. 27 
in Eq. [26] yields 

R- 0 3 \ / pdi \ , 


Se di = M (ex) 


0 , 


V* 


1 0 0 0 0 0 Se di = 0 V i = 1 . . . i 


(28) 


where R is the rotation matrix defined in Eq. |27[ 0 3 is a 3 x 3 matrix of zeros, and I 3 is a 3 x 3 identity 
matrix. The variables being chosen in the new optimization are the vectors 5e d i, the rotations 0 X , 0 y , and 
0 Z (contained in R), the velocities v d i, and the translations t. The optimization can perform small rotations 
(0 X , 0 y , and 0 Z are constrained), but can no longer optimize the tetrahedron scale and still retain linearity. 


Example: Fuel Expenditure for Rotation Approach — Using the highly elliptical orbit from the 
example in Section |III[ the initial conditions for a tetrahedron with 1000 kilometer sides were found using 
the optimization criterion in Eq. [25] and the constraints in Eq. [28] and Eq. [26] It was assumed that the four 
satellites began from the near-tetrahedron configuration 
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Table 1: Maneuver Cost: Tetrahedron Initial Condition Optimization 


DOF 

SC 1 

SC 2 

SC 3 

SC 4 

Initialization Total 

No Rotation 
No Translation 

51-fl 

193 

213 

6.487 x 10 6 

6.487 x 10 6 

No Rotation 
Translation: ±1 km 

127 

246 

161 

299 

833 

No Rotation 

Scaling: ±5%, Translation: ±1 km 

101 

423 

0.430 

202 

726 

Rotation: ±0.1 radians 
No Translation 

8.316 

7.696 

5.846 

1.484 

23.342 

Rotation: ±0.1 radians 
Translation: ±1 km 

2.422 

0.669 

1.080 

1.199 

5.370 


When computing the objective function for this example, the spacecraft state errors and individual elements 
were weighted equally. Table [T] shows the results of several optimizations, each with different degrees of free- 
dom enabled. The table shows that the most significant fuel advantage is achieved through the combination 
of rotation and translation. In the rotation/translation combination, the optimization created the desired 
differential elements 


(30) 


where the limits of the translation variables are ± 1 km and the limits of the three rotation variables were 
each ±0.1 radians. Within those limits, the optimal translation and rotation vectors were found to be 

t = ( 26.809 52.500 -6.344 ) T (31) 

( d x d v 9 Z ) T = ( -0.001005 -0.000997 -0.0007217 ) T (32) 
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in meters and radians. The planner described in Section IV was used to create the formation geometries 


specified in Eq. 30 from the initial conditions specified in Eq. 29 The total maneuver fuel cost to create the 
formation ( i.e the sum of the costs for each spacecraft) was 6.33 mm/s. 


VII. Formation Flying: Coordination Using GVEs 


The model predictive controller described in Section |IV| is designed to be decentralized, with a fully 
independent controller being run on each spacecraft. The controller designs trajectories that will keep 
a spacecraft i inside an error box centered about the spacecraft’s desired orbit, e*. In Section III the 
desired orbits are defined with respect to the actual orbit of an arbitrary satellite in the formation, ei, 
using differential orbital element vectors, Se^, in the same manner used in Ref. [13j Section |VT| presented 
an approach for choosing initial conditions that minimized the weighted state error across the formation. 
In a system where initial conditions are chosen infrequently, it may be desirable to introduce additional 
coordination into the formation. When spacecraft each track desired states with no coordination, the control 
task is referred to as formation-keeping . 5 Alternately, formation flying occurs when the spacecraft controllers 
collaborate to achieve formation-wide fuel minimization. This coordination can be achieved by calculating a 
central point that minimizes the overall weighted state error of each spacecraft in the formation. Approaches 
to implementing closed- loop coordination of this type are presented in Refs. 0 and |22[ The virtual center 
approach in Ref. [Tl] is a centralized calculation of the error-minimizing center based on fuel-weighting and 
derived from measurements available through carrier-phase differential GPS (CDGPS) relative navigation 


“Fuel use indicated in mm/s 
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of the type described in Ref. 23 An equivalent approach can be used to find an error-minimizing reference 
orbit for a formation described in differential orbital elements. 

Measurements from a CDGPS relative navigation system are assumed to be in the form of relative LVLH 
states [23} x, : , for each satellite in the formation. The measurements will be relative to an arbitrary absolute 
satellite state, ei, in the formation, which is assumed to be at the origin of the LVLH frame. In addition 
to relative states, the GPS sensors on each satellite can be expected to compute a less accurate estimate of 
the spacecraft’s absolute state. Given an estimate of the absolute state in Earth Centered Inertial (ECI) 
coordinates, Xecp and the relative states x M the differential states 5e in Eq. [4] can be computed in several 
ways. The matrix M(e i) from Eq. 22 could be computed and used to create a first order approximation of 
the relative differential element states. However, an exact conversion can be calculated by forming estimates 
of the absolute states of each of the satellites based on their relative measurements 


X^eci; = Xecix + Xj (33) 

The absolute states Xecp can be converted to Keplerian orbital elements, e,;, of each satellite using a well- 
known procedure described in Ref. [l] The relative measurements are then recovered in terms of differential 
orbital elements, <5ej, using Eq. [4j Desired differential elements, Se dd, are then specified with respect to 
an unknown virtual center state, (5e c , which is specified with respect to the absolute state ei. Using the 
procedure described in Ref. 11 the error of a spacecraft with respect to the virtual center, (d is given by 

e i ()O r (d (31 ) 


which can be placed in the standard least squares form bi — A^e c = (d, where bi = e 7 — Sedd, Ai is a 
6x6 identity matrix, and 5e c denotes the location of the virtual center with respect to ej in differential 
orbital elements. By concatenating the Aj, and ( ci vectors for each spacecraft, the statement of error 
for the entire formation is written b — A8e c = (, where b = ( b\ ... b n ) T , A = ( A\ ... A n ) T , and 
( = (( c i . . . (cn ) T . The solution that minimizes the error vectors globally in a weighted least squares sense 

1S <5e c = ( A T WA)~ 1 A T Wb (35) 


where W is a weighting matrix that can be used to bias the center location according to the fuel-use rates of 
different satellites in the formation, as well as to weight orbital elements individually based upon the amount 
of control required to alter them (obtainable from the GVEs for ei). 


VIII. Decentralization of Virtual Center Scheme 

The virtual center calculation in Eq. [35] is the well-known closed- form solution to a least squares problem. 
The solution to the decentralized least-squares problem is also readily available in the literature and is given 
for spacecraft i as 

Se Ci = 5e Ci _ 1 + (A T W i A)~ 1 AjW i (b i - AiSe Ci _ J V i = l...n (36) 

where A = ( Ai ... Ai ) T and IV, = diag (Wi, . . . , W*). This decentralized procedure would begin with 
an initial estimate <5e Cl = bi and spacecraft i would pass both W, and 8e c . on to the next spacecraft (i + 1) 
in the formation. The estimate of the virtual center, Se n will then be equal to the fully centralized solution, 
which could be broadcast back to the entire formation. 

The decentralized calculation can be simplified using several reasonable assumptions. First, assume that 
the relative weighting between differential orbital elements is the same for all spacecraft and is calculated 
to be W e based on the absolute orbit which is known by all spacecraft in the formation. Next, let the 
fuel-based weighting for each spacecraft be represented by a scalar, w^. Now, the total weighting matrix, 
W, , for each spacecraft is given by 

Wi=WiW e V i = l...n (37) 

Also assume that the variable being solved for is always the exact virtual center, Se c , making the A* 
matrix for each satellite a 6 x 6 identity matrix. Now, for the ith satellite W = diag (w\W e , . . . ,WiW e ), 
and A is a column of i concatenated 6x6 identity matrices. The decentralized solution for the itli satellite 


10 of H2I 


2nd International Symposium on Formation Flying Missions & Technologies 


becomes 


Se Ci = <5e Ci _ 1 + (A t WA) 1 AjW i {b l -A i 5e Ci _ 1 ) 

( 


= Se Cil +Wi 


I ... I 


w\ W e 0 ... 0 

0 vj‘2 W e ... 0 
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WiW R 
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i 
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W e (bi - <5e Ci _J 


which can be further simplified to 

Se Ci = Se c ._ 1 + w i [{w 1 + . . . + w i )W e ]~ 1 W e {b i - <5e c ._J = Se Ci _ 1 H (6» - (5e Ci _J 

W\ + . . . + Wi 

Introducing u>i, the sum of all previous scalar weights Wi = J2k=i w k, the new recursion becomes 

~ - ( bi 

Wi - 1 + m 


(38) 


(39) 


(40) 


and the only information that needs to be passed to the next spacecraft to form <5e Ci+1 is the current virtual 
center estimate, 5e Ci , and the scalar u>i. 


IX. Formation Maintenance on MMS-like Mission 


The control system described in Section |IV| was 
demonstrated on a segment of the MMS mission. The 
MMS mission is comprised of four spacecraft that cre- 
ate regular tetrahedron geometries once per orbit. The 
orbits of the four spacecraft are widely separated and 
highly elliptical, presenting a challenge for many opti- 
mal formation specification and control approaches in 
the literature 00 Using the tetrahedron initial-condition 
optimization approach in Section VI and the model pre- 
dictive approach in Section [IVj the four spacecraft were 
controlled in a fully nonlinear simulation with atmo- 
spheric drag effects (unmodeled oblateness disturbances 
were not used) using a commercial orbit propagator. 
Figure [3] shows the rate at which fuel was used over the 
course of two weeks of formation flying. After an initial 
transient period, the formation fuel use rate converges 
to approximately 5.95 mm/s per day (« 1 orbit) for each 
satellite. 


Formation Flying in a 1000 km Tetrahedron 



Fig. 3 : Fuel cost for maintaining a 1000 km tetra- 
hedron formation in a highly eccentric orbit 


X. Conclusion 

Gauss’ Variational Equations have been used to derive a set of linearized relative dynamics of orbital 
motion. Using a linear parameter varying dynamics model, such as GVEs, for control allows a compromise 
between simple but inaccurate linear models (e.g., Hill’s equations) and high fidelity but difficult to control 
nonlinear models. The linearization assumptions for the new set of dynamics were shown to be valid for 
typical spacecraft error box sizes. The linearized GVE-based dynamics were used in a model predictive 
controller of the form described in Ref. [5] and the combination is shown to be more fuel-efficient than an 
impulsive based (non-optimized) GVE technique. By combining differential element orbit specifications and 
dynamics linearized about the spacecraft in the formation, the problems of long-baseline initialization and 
control discussed in Ref. [24] are eliminated. 

A method of specifying differential orbital elements for a desired formation geometry was introduced which 
optimizes the parameters in the initialization problem to minimize the fuel required to maneuver into the new 
formation. This method was applied to a regular tetrahedron-shaped orbit of the type planned for the MMS 
mission. A method of applying rectilinear error box constraints to a formation specified in differential orbital 
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elements was also presented. Finally, the virtual center approach for coordinating spacecraft format iond^ 
was formulated in terms of the differential orbital elements, and a means of decentralizing the algorithm was 
shown. The GVE-based dynamics/MPC controller was used to specify and control a tetrahedron-shaped 
formation in an MMS-like orbit for a period of two weeks. 
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